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|]  ABSTRACT 

An  efficient  contour-plotting  routine  is  discussed  which 
is  based  on  a  scanning  algorithm  of  Cottafava  and  LeMoli  and 
employs  bi-linear  interpolation.  An  auto-interpolation  scheme 
is. developed  which  automatically  adjusts  the  number  of  inter¬ 
polations  in  any  data  square  to  produce  smooth  line  segments. 
A  program  listing  and  examples  are  given. 
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ABSTRACT 


An  efficient  contour-plotting  routine  is  discussed  which 
is  based  on  a  scanning  algorithm  of  Cottafava  and  LeMoli  and 
employs  bi-linear  interpolation.  An  auto-interpolation  scheme 
is  developed  which  automatically  adjusts  the  number  of  inter¬ 
polations  in  any  data  square  to  produce  smooth  line  segments. 
A  program  listing  and  examples  are  given. 
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INTRODUCTION 


There  are  many  approaches  to  producing  contour  maps  on  a 
digital  computer.  Several  of  these  are  described  by  Cotta,aVa 
and  LeMoli  (1969).  Ideally,  a  contour  routine  should  be  e  fi  "t 
on  both  computer  ahd  plotter,  but  most  seem  to  possess  only  on . 
kind  of  efficiency.  For  quite  a  number  of  reasons,  inc  u  ing 
plotting  efficiency,  algorithms  of  the  -line-following  type 
are  preferable  for  use  with  mechanical  plotters.  Cot  a  ava 
LeMoli  present  a  scanning  algorithm  of  this  type  which  is  «  so 
very  efficient  on  the  computer. 

in  their  program  (Cottafava  and  LeMoli.  private  communica¬ 
tion)  thoy  assumed  a  linear  variation  between  points,  but  did 
not  do  any  interpolation  in  the  interior  of  the  data  square 
(defined  by  four  contiguous  data  points  as  vertices).  Thus  a 
plot  produced  with  their  routine  consists  entirely  of  straight 
line  segments  joined  together,  the  coarseness  depending  on  the 
spacing  of  data  points.  To  remedy  this,  the  author  *••!>» •" 
autointerpolation  scheme  employing  bi-linear  interpolation 
use  in  the  Interior  of  the  data  square.  With  this  scheme, 
described  in  the  present  paper,  the  number  of  Interpolations 
used  in  crossing  the  data  square  is  automatically  adjusted  to 

produce  a  smooth  curve,  the  number  required  depending  on  th* 
curvature  of  the  line  segment.  This  method  requires  no  additional 

storage  and  is  very  fast. 

While  the  interpolation  does  indeed  produce  smooth  line 
segments,  slope  discontinuities  sometimes  occur  on  the  edges 
(of  the  data  squares).  This  is  a  limitation  of  the  bi-lin 
interpolation  law.  An  easy  and  economical  solution  to  this 
problem  is  to  obtain  a  finer  data  mesh  by  performing  a  h,g  e  - 
order  interpolation  before  entering  the  contour  routine.  This 


may  also  be  done  when  the  data  are  not  given  as  a  set  of 
equally-spaced  points. 
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SCANNING  ALGORITHM 


The  scanning  algorithm  used  is  that  of  Cottafava  and 

wlth  “''"i"  m0£f'ificat1ons  by  this  author.  A  more 
complete  discussion  than  will  be  given  here  can  he  found  in 
the  paper  cited.  For  each  contour  value,  the  procedure  Is  to 
scan  the  entire  data  array  to  find  which  line  segments  are 
intersected  by  the  level  line  and  to  store  the  Information  as 
flags  within  the  data  words  themselves.  A  horizontal  and 
vertical  segment  are  associated  with  each  data  point  as  shown 

'9U,"e  1  (the  J'-axis  Siven  its  normal  sense  here- 
Cottafava  and  LeMolt  reverse  it). 


^ _ A. 


Figure  1.  Intersection  flags 


The  flags  are  stored  as  bits  in  the  upper  part  of  the  Integer 
word,  the  data  being  normalized  to  positive  Integers  by  a  Hn 


transformation.  (In  deference  to  readers  with  machines  which 
cannot  perform  masks,  the  use  of  masks  in  the  programs  has 
been  limited  to  an  inessential  role,  where  they  can  easily  be 
replaced  by  a  test  and  subtract.  Where  masking  statements  are 
available,  the  flags  can  be  conveniently  stored  as  the  least- 
sigriificant  bits  of  the  floating-point  mantissa  of  each  data 
word.  This  results  in  a  considerable  simplification  of  the 
program. ) 

The  contours  are  traced  in  a  second  scanning  operation, 
each  flag  being  erased  as  it  is  found.  This  procedure  makes  it 
easy  to  find  all  the  branches  of  the  contour  level  and  is  in 
large  part  responsible  for  the  efficiency  of  the  program.  Each 
line  is  traced  square  by  square  by  a  local  scan  which  checks 
all  the  edges  of  the  data  square  in  fixed  order  (counter¬ 
clockwise  beginning  with  the  right  edge)  to  find  the  contin¬ 
uation  of  the  line.  This  procedure  runs  into  trouble  only  in 
the  case  of  an  interior  saddle  point  (square  crossed  twice  by 
the  same  contour),  and  in  this  case  there  is  an  easy  solution 
based  on  the  interpolation  method. 
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INTERPOLATION 


In  order  to  define  the  behavior  of  the  contour  line 
inside  each  data  square,  we  must  make  an  assumption  about  the 
behavior  of  the  function  inside  the  square.  Lacking  any  special 
information  in  the  general  case,  we  assume  bi-linear  variation 
as  the  simplest  general  variation.  (Using  a  higher  order  inter¬ 
polation  here  would  also  cause  serious  difficulties  with  the 
scanning  procedure.)  That  is,  we  assume 


F(x,y)  =  A  +  0x  +  ay  +  <5xy  (1) 

referred  to  a  local  coordinate  system  with  origin  at  A,  as 
shown  in  Figure  2. 

Y 


t  Y 

taj tc  > - - o 


Figure  2.  Interpolation  conventions 


Equation  (1)  is  the  Lagrange  2x2  interpolation  formula  and 
is  also  equivalent  to  a  Taylor's  series  expansion,  with  the 
assumption  of  constant  first  derivatives  on  each  edge.  Making 
our  square  of  unit  dimensions,  the  derivatives  are 


a  =  C-A 
3  =  B-A 
y  2  D-C 

6  2  y-e 


With  these  definitions,  it  is  easy  to  see  that  (1)  reduces  to 
the  correct  values  at  the  corners  and  reduces  to  ordinary 
linear  interpolation  on  each  edge.  Our  contour  segment  is 
therefore  the  locus  F(x,y)  -  v,  the  value  of  the  contour. 

From  (1)  we  obtain  either 


y 


v-A-Bx 

a+6x 


(3) 


or 


_  v  -  A  -  ay 
3+<5y 


(4) 


with  0<x<1  and  0<y£l .  These  expressions  are  easily  computed. 

In  practice,  we  choose  the  number  of  interpolations,  subdivide 
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Ax  or  Ay,  and  use  either  (3)  or  (4),  respectively.  It  is 
convenient  to  use  (3)  for  a  line  terminating  on  a  vertical 
edge,  and  (4)  for  a  line  terminating  on  a  horizontal  edge. 


AUTO-INTERPOLATION 


We  approximate  our  ideally  smooth  curve  by  a  series  of 
chords.  If  we  choose  the  number  of  chords  in  each  square  so 
that  the  maximum  deviation  from  the  ideal  curve  is  on  the 
order  of  the  basic  plotter  increment,  then  we  obtain  as  smooth 
a  curve  as  we  can  with  no  wasted  time.  We  take  as  our  auto¬ 
interpolation  criterion,  the  maximum  perpendicular  distance 
from  the  curve  to  the  straight  line  between  the  end  points. 

The  perpendicular  distance  from  a  point  to  a  line  is 

e ( x )  =  (y  -  mx-b)/Ym2+l  (5) 

from  elementary  geometry.  We  consider  the  case  of  a  segment 
terminating  on  the  left  edge  of  the  square,  and  for  this  case 
there  are  just  two  distinct  possibilities,  as  shown  in  Figures 
3  and  4.  All  other  possibilities  can  be  found  by  reflections 
and  a  rotation. 


c  &  *  a  c  o 


A  S  A  B 

Figure  3.  Case  I  Figure  4.  Case  II 
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For  case  I  the  slope  is 


v-C  ,  v-C  _  , 

nil  _  C-A  1  C-D  "  ~y/a 


For  case  II, 


v-A-3  _  v-A  _  _  6(v-A)  +  a3 
m2  a+6  a  a ( a+6 ) 


For  both  cases  the  intercept  is  b  =  (v-A)/a.  The  "brute  force" 
calculation  of  the  maximum  of  (5)  is  messy,  hut  a  transformation 
simplifies  the  algebra.  Define 


n2  =  6/a 

n  h  —  Ax  where  Ax 

a 


(v-C)/Y,  case  I 
1 ,  case  II 


(6) 


Then 

+  B  =  6[1  +  +  s  ■  y[l  +  n,] 
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which  leads  to  the  relation 


m2  (l+n2)  =  mi  (1+r>] ) 


It  can  also  be  shown  that 


miO+n-i) 

y_5  =  - ‘ —  x 

1  +  n  2  x 


and 


(7) 


y' 


™  (l+n2x)2 


The  location  of  the  maximum  of  (5)  is  given  by  the  condition 


e'(x)  =  0  =  y ' ( x)  -  m 

(hats  will  be  used  to  refer  to  the  maximum),  i.e., 
m-|  (l+n1 ) 

- —  =  m 

( 1 +n2  x ) 2 


(8) 


-10- 


Thus 


y-b  =  mO+r^*)* 


and 

V-2  ,  a  a  ,  a  a2 

.  rrr+1  e=y-b-mx=  mn  g  x 


or 


A 

E 


;2 


where  x  is  found  from  (8).  For  case  I,  m  =  m-j ,  and 


(9) 


*1  ■  k  ^ 

For  case  II,  m  =  m2>  and  using  (7)  in  (8)  gives 
*2  =  7^  ‘  '] 


which  has  the  same  form.  Putting  these  results  in  (9)  gives 
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for  both  cases.  This  is  conveniently  written  in  the  form 


<S«/rer-iM  no) 

VAx  +  Ay 

Since  (10)  is  symmetric  in  x  and  y,  it  can  be  easily  seen  to 
apply  to  the  cise  of  a  segment  terminating  on  the  bottom  edge 
of  the  square  also,  provided  that  we  define  n  =  6Ay/3  in  this 
case.  Thus  all  possible  cases  are  contained  in  (10).  A  bound 
can  be  placed  on  e  by  considering  it*00  and  Ax-*Ay->-l ,  namely 
e<l //~T,  which  agrees  with  geometrical  intuition. 

Knowing  how  to  calculate  the  maximum  deviation  of  the  curve 
from  the  straight  line  between  endpoints,  we  use  this  to  estimate 
the  number  of  segments  needed  to  approximate  the  curve  to  any 
desired  accuracy.  To  do  this,  we  consider  the  case  of  a  circular 
arc.  Figure  5. 


Figure  5.  Deviation  from  a  chord 


It  is  easily  shown  that  dsL  /8r,  provided  L/r<<l,  the  important 
point  being  the  proportionality  d  «  L^.  Since  L  “  1/N,  approximately, 
where  N  is  the  number  of  segments  needed  In  that  square,  we  take 


N  «  1  +  /]T]7I 


where  d  equals  the  allowable  deviation.  For  plotters  with  a  2.5  mil 
increment,  d  =  .001"  is  satisfactory.  This  auto-interpolation 
scheme  appears  to  work  quite  well. 


SINGULARITIES  AND  SADDLE-POINTS 


Evidently  something  peculiar  happens  with  (10)  if  n+l<0. 
This  can  be  seen  to  be  connected  with  a  singularity  in  y  or  x, 
equation  (3)  or  (4).  The  existence  of  a  singularity  in  y(x) 
within  the  interval  0£x£l  is  implied  by  the  condition  that  (D-B) 
and  a  have  opposite  signs,  and  the  existence  of  a  singularity  in 
x(y)  within  the  interval  0_<y£l  is  implied  by  the  condition  that 
Y  and  3  have  opposite  signs.  If  both  singularities  exist,  then 
the  square  has  an  interior  saddle  point  located  at  the  inter¬ 
section.  An  example  is  shown  in  Figure  6.  A  very  conveni ent 
criterion  for  making  connections  in  a  saddle  square  is  that 
contours  should  never  cross  a  singulari tv.  The  condition 
n+l<0  can  be  easily  shown  to  imply  a  wrong  connnection  in  a 
saddle  square.  By  far  the  easiest  solution  to  this  problem  is 
just  to  check  for  a  wrong  connection  and  to  resume  scanning 
the  square  if  it  exists.  At  most  three  tries  will  be  necessary 
to  make  the  correct  connection,  and  since  saddle-points  should 
be  relatively  rare,  this  is  a  small  price  to  pay  for  such  a 
simple  procedure  that  guarantees  correct  connections. 


*o  =-a/S=4 

-  /3/S  -  J'g 


Figure  6.  Contours  in  a  saddle  square 
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RECTANGULAR  MESH 


Should  it  be  desired  to  plot  data  calls  as  rectangles 
instead  of  squares,  this  is  easily  done  by  stretching  one  axis 
in  the  calls  to  PLOT.  Defining  r  as  the  ratio  of  x  to  y  scale 
factors,  i.e.  the  rectangle  has  length  r  in  the  x-direction 
and  1  in  the  y-direction,  elementary  trigonometry  gives  for  the 
modified  deviation 


e  =  e 


-,/^V 

\l^~ 


\ 

i 

t 

! 
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EXAMPLES 


Two  simple  examples  of  plots  produced  by  the  routine 
are  shown  in  Figures  7  and  8.  Figure  7  was  produced  from  real, 
deterministic,  data,  and  Figure  8  from  random  numbers.  Data 
points  are  marked  by  ticks  along  the  borders  of  the  plots.  Each 
plot  is  based  on  only  24  data  points,  and  the  large  number  of 
slope  discontinuities  indicates  the  need  for  a  more  refined 
data  mesh. 
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Figure  7.  Example  1,  real  data 


Figure  8.  Example  2,  random  numbers 


APPENDIX 


PROGRAM  LISTING 

Written  in  FORTRAN-63,  a  programming  language  of  the 
CDC  1604  computer.  Integer  words  assumed  at  least  33  bits  long. 
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IS  ■  I  ■  I  -  1 

IF  <  I v  < I # J*i )  .Ug,  I  VfcR  )  10*15 

4  a  J*x  %  H|SP  ■  2  *  ASS  I  UN  5  10  JAIL  S  QO  TO  20U 

IbV  ■  mi*t.J».ANO..NOT.lVfcH 
IF (  IbV.Rfc.lHOH  >  20.25 

|  a  |*i  f  M I SP  ■  1  A  ASSIUN  5  >0  JAIL  A  QO  To  200 

IF (  I  V 1 1 # J > ,Gfc* I VfcR  >  3q.3» 

N 1 SP  >2  A  ASSIGN  3  TO  JAIL  A  QO  TO  20U 

IF (  I  V ( I , J ) , Gt »  I H0R  )  4Q.90 

N ISP  *  i  *  ASSIGN  4  TU  JAIL  *  QO  TO  2qu 

IF(  K  >  ioo»l90 

IlfcTURN 


fcMOVfc  flags.  iNTfcHPOLATP.  ANO  PLOI  POINT 
.  ihfcrl  A  J  ■  jfigM  A  IVCI.JJ  »  IVU#J»>1VEH 

ISP  ■  2  *  ASSIGN  mg  TO  JAIL 

0  ■  0  •  2  *  VO  ■  I  •  2 

0  i  dP  A  10  ■  IP 


FI  J.fcQ.JP  >  210.220 
0  ■  JP  ♦  1 

FI  I. fcO. IP  )  23o»2«u 
0  ■  IP  ♦  1 

■  IA  ■  lVII*J).ANO..NOT(iVH 

a  lb  •  ivii»jq)*and..not*ivh 

a  It  a  IVIIU»J).AND..NOTiIVH 
a  It  a  IVtlU*JQ».ANO..NO!«IVH 
LF  ■  0  -  A  A  UET  a  B  ”  a  A  DbL  a  U-C  *  BfcT 

0  TO  l250*2Nb^»  NlSP 

F  a  IV0-A)/PET  ♦  XO  A  0*  «  *F  •  IP  *  YF  « 

IIGN  a  io  .  i  A  OT  ■  V  a  SIUNtl VP9TO) 

TA  a  OY*OEL  •  BET  »  QO  10  270 

F  ■  IVO-AI/aLF  ♦  YO  A  DT  •  YF  -  TP  $  XF  ' 

IGN  •  JO  -  J  *  0*  ■  *  ■  SIGN*<XP-X0) 

TA  a  UX*DEL  i  AfcF 
FI  ETA  )  28n»400 
TAP  ■  ETA  ♦  i. 

ADDLfe-POlNT  TfcST 
Ft  ETAp  )  290*295.300 

I  |S  A  J»JS  *  IT"IT4| 

0  TO  I .5*25.35 ) •  IT 

PS  a  1.  A  GO  To  iln 

PS  •  DX*DY*ISURTM6TAP)»1*I**2  /  tfc  \  A*SUHTF  tUX*UX*DYtDY  >  > 

I  •  1.5  ♦  SQHTFIABSFIEPSI'UEVS) 

NI  a  N|  *,99W  A  IT  •  || 

0  TO  t3?0*35i|)*  TIISP 


» 

» 


UP  ■  -ALP*J* 
Jb  ■  UtL*UY 


32b 

330 

840 

35o 

B60 

670 

400 


QY  ■  Y  /  F  N 1 
*  a  VO-A  -  ALF*Y 
8  a  Btt  ♦  DEL** 
t  a  y  -  nr 

IF  (  Y  )  4U0#4'tU,j4o 
S  a  f  •  flf  *  <i 


UG 

1P6N 


OALL  BLOT <  XO*X* YO*SIGN*Y j 
QO  TO  330 

IX  ■  x  /  fni 

X  a  VC-A  -  «pT*X  * 

Q  a  ALP  ♦  DFL«X  *  D1* 

1  a  x  -  nx 

IF  C  X  >  40b#4JU#370 

P  a  F  •  Of  *  G  ■  6  •  UG 

8ALL  PLOT  <  X0*S1GN*X. YO*Y#  IP£N  1 

QO  TO  36 C 

UMI.jI  ■  -  1HORABI&P 

IS  *  IP  ■  I  *  JS  »  JP  *  J  5 
*ALL  PLOT  <  XP#  TP*  IPEN  ) 

QO  TO  JAIL*  <X*3»4#5*100> 

IND 


X  ■  Mb 


“BE  I  *UX 
UEL*UX 


F/G 


XP  ■  XP 


